Finite temperature phase diagram of spin- 1/2 bosons in two-dimensional optical lattice 
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We study a two-species bosonic Hubbard model on a two-dimensional square lattice by means of 
quantum Monte Carlo simulations and focus on finite temperature effects. We show in two different 
cases, ferro- and antiferromagnetic spin-spin interactions, that the phase diagram is composed of 
a superfluid phase and an unordered phase that can be separated into weakly compressible Mott 
insulators regions and compressible Bose liquid regions. The superfluid-liquid transitions are of 
the Berezinsky-Kosterlitz-Thouless type whereas the insulator-liquid passages are crossovers. We 
analyse the pseudo-spin correlations that are present in the different phases, focusing particularly 
on the existence of a polarization in this system. 
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I. INTRODUCTION 

The study of strongly interacting quantum models by 
direct realization of an experimental system reproduc- 
ing the model properties, an idea proposed by Feynman 
[H, was realized in the past ten years with the produc- 
tion of Bosc-Einstcin condensates (BEC) and their use 
as "quantum simulators" Q. In particular, BEC in con- 
junction with optical lattices are used to produce systems 
reproducing the physics of well known quantum statisti- 
cal discrete models such as fermionic or bosonic Hubbard 
models. 

Used to study simple models of bosons [j| or fermions 
Q at low temperature, the flexibility offered by these 
systems extends the range of interesting models to more 
exotic ones which can be treated both experimentally 
and theoretically. Examples include systems with long 
range interactions Q , fermions with imbalanced popula- 
tions Q, mixtures of different kinds of particles and 
spin-1 bosons with spin- independent and spin-dependent 
interactions which allow interplay between superfluid- 
ity and magnetism [8MTol| . Furthermore, it is possible 
to study systems of bosons with two effective internal 
degrees of freedom on an optical lattice, the so-called 
"spin-1/2 bosons". Such a system, with spin-dependent 
interactions, could be produced by applying a periodic 
optical lattice on a bosonic system with two triply de- 
generate internal energy levels. The optical potential ap- 
plied would localise the atoms at the nodes of a regular 
network, but would also couple the internal states by A 
or V virtual processes, thus leaving only two internal low 
energy degenerated states denoted and A and realiz- 
ing an effective spin-1/2 model [HI, The presence of 
the spin-dependent interaction introduces a term in the 
Hamiltonian which permits the conversion of two parti- 
cles of one type into the other type and renders numeri- 
cal simulations more difficult. This model is related, but 
not identical, to other models including p-band super- 
fluid models [l3l - fl5j and the bosonic Kondo model fltjl ]. 
Understanding the phase diagram and properties of the 



simpler spin-1/2 bosonic system takes us a step toward 
understanding the more elaborate, and more difficult to 
simulate, models. 

The spin-1/2 model has been extensively studied with 
mean-field theory (MFT) at zero or finite temperatures 
and, in previous work, we explored its zero temperature 
behavior with quantum Monte Carlo (QMC) simulations 
in one [l?] ] and two [l8| dimensions for on-site repulsive 
interactions. At zero temperature, the phase diagrams 
obtained in one and two dimensions, with MFT or QMC, 
are similar. Generally speaking, at zero temperature, the 
system can adopt two different kinds of phases: insulat- 
ing Mott phases that appear for integer density, p, and 
for large enough repulsion between particles and super- 
fluid phases (SF) otherwise. The detailed nature of these 
phases depends on the interactions between the differ- 
ent kinds of particles. In the case where the repulsion 
between identical particles is smaller than between dif- 
ferent particles (U2 > in our previous work fl8| and in 
the following) the superfluid is found to be polarized, that 
is, an imbalance develops in the populations of the two 
kinds of particles and one of the species becomes dom- 
inant. The p — 1 Mott phase is also polarized whereas 
the p = 2 phase is not (we did not study higher densities 
with QMC). In the opposite case, U2 < 0, all the phases 
are unpolarized. Noteworthy is the presence of coherent 
exchange movements ffiSj , where two particles of different 
types exchange their position in the p — 1 Mott phase 
in both cases, as well as in the p = 2 Mott phase for 
U2 < 0. Finally, in one dimension, all zero temperature 
phase transitions were found to be continuous whereas 
in two dimensions and when U2 > is small enough, the 
p = 2 Mott-superfluid transition was predicted by MFT 
to be first order near the tip of the Mott lobe and con- 
tinuous otherwise, whereas for larger U2 the transition 
was predicted to be always continuous. This was con- 
firmed by QMC simulations Related spin-1 models 
were studied using MFT [2fj|. l2f| or QMC in one dimen- 
sion [Io| and a similar spin 1/2 bosons model was also 
recently studied 22 1. 
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In this paper, we will study the spin-1/2 model at finite 
temperature in two dimensions and compare with MFT 
predictions. The results for finite system sizes and tem- 
peratures are relevant to experimental efforts to study 
this system. In Section II, we will introduce the model 
and the MFT and QMC techniques used to study it. Sec- 
tion III and IV will be devoted to the presentation of the 
results obtained for the U2 > and U% < cases, respec- 
tively. We will summarize these results and give some 
final remarks in Section IV. 



II. SPIN-1/2 MODEL 



The model we will study is the same we previously 
studied in the low temperature limit in one 1171 and two 
dimensions f]~sj | and previously introduced in |ll| . It is an 
extended Hubbard model governed by the Hamiltonian 
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"T E ( a 0r a 0r a ArOAr + a Ar a Ar a 0ra0r ) , (3) 



where operator a ar {a aT ) destroys (creates) a boson of 
type a = 0, A on site r of a two-dimensional square lattice 
of size L x L. The operators h ar measures the number of 
particles of type a on site r. The densities of particles of 
type and A are called po and p\ while the total density 
is called p = po + pa- 

The first term |T]) of the Hamiltonian is the kinetic 
term that lets particles hop from site r to its nearest 
neighbours r'. The associated hopping energy t = 1 sets 
the energy scale. A chemical potential p is added if one 
works in the grand canonical ensemble. The second term 
^ describes on-site repulsion between identical particles 
with a strength Uq or between different particles with a 
strength Uo + l^- We will study both the positive and 
negative U2 cases but will keep only repulsive interac- 
tions, that is |Z7a| < Uo, and a fixed moderate value of 
\U2\/Uq = 0.1. The last term ([3]) provides a possibility to 
change the "spins" of the particles: When two identical 
particles are on the same site, they can be transformed 
into two particles of the other type. It was demonstrated 
that the matrix element associated with this con- 
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version is U2/2. We are using a different sign for the term 
([3]) compared to the articles where the model was orig- 
inally introduced fll| but we have shown in a previous 
work (ItI ] that this sign can indeed be chosen freely due 
to a symmetry of the model. 



A. Mean Field Theory 

The only term that couples different sites in the Hamil- 
tonian is the hopping term ([T]). Introducing the field 
ip a = (a, aT i) = (a<Tr'}, we replace the creation/destruction 
operators on site r' b y t heir mean values ip a , following 
the approach used in [llj. The Hamiltonian on site r is 
then decoupled from neighboring sites and can be easily 
diagonalized numerically in a finite basis. The optimal 
value of the fields ip a are then chosen by minimizing the 
grand canonical potential Q = — feTln Q with respect to 
ip a where Q is the grand canonical partition function. 
The system is in a superfluid phase when ip c is nonzero, 
with superfluid density p s = \ipo\ 2 + |"0a| 2 , and is other- 
wise in an unordered phase where two cases can be dis- 
tinguished: an almost incompressible case, i.e. a Mott 
insulator, and a compressible case, i.e. a liquid. 

In these two cases there is no broken symmetry and 
they cannot be distinguished by symmetry considera- 
tions. If there is a first order transition between the 
MI and the liquid, characterized by discontinuities in 
the density or other thermodynamic functions, the MI 
and the liquid would be two distinct phases. If, how- 
ever, there is no discontinuity in the evolution from the 
MI to the liquid then they are only two limiting cases 
of the same unordered phase and there is only a smooth 
crossover between the MI and the liquid regions of the 
phase diagram. We shall see that, indeed, there is a 
crossover in the system we are considering here. 

One can distinguish between almost-incompressible 
and liquid regions by calculating the local density vari- 
ance which is a measure of the local compressibility, 
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(n r ) ) where n r is the total number of 



particles on site r. k is close to zero in the Mott phase 
and much larger in the liquid phase. 

While this MFT was shown to reproduce qualitatively 
the phase diagram at zero temperature [l8[, it is rather 
limited at finite temperature. Indeed, whenever the ip G 
are zero the hopping parameter t no longer plays a role 
in the MFT. Then, while the MFT can distinguish be- 
tween SF and unordered phases, it does not correctly dis- 
tinguish Mott Insulator (MI) regions from normal Bose 
liquids ones, as the crossover boundary between those re- 
gions will not depend on t and will be the same as in the 
t = case. 



B. Quantum Monte Carlo simulations 

To simulate this system, we used the stochastic Green 
function algorithm (SGF) [23|, an exact Quantum Monte 
Carlo (QMC) technique that allows canonical or grand 
canonical simulations of the system at finite tempera- 
tures as well as measurements of many-particle Green 
functions. In particular, this algorithm can simulate ef- 
ficiently the spin-flip term in the Hamiltonian. We stud- 
ied sizes up to L = 14. The density p is conserved in 
canonical simulations, but individual densities po and p\ 
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fluctuate due to the conversion term Eq. (|3]). The super- 
fluid density is given by fluctuations of the total winding 
number, (Wo + W\), of the world lines of the particles 


((Wq + W A ) 2 ) 

Ps = — Wp — ■ (4) 

The superfluid density cannot be measured separately 
for and A particles due to the conversion term [2]|. It, 
therefore depends on the total density not on the separate 
densities of the two species. We also calculate the one 
particle Green functions 

Ga ^ = 21? ^ ( a lr+R a *r + a lr a *r+R) > ( 5 ) 
r 

with a = 0,A. G CT (R) measures the phase coherence of 
individual particles. 

In a strongly correlated system it is useful to study 
correlated movements of particles which can be done, 
for example, by studying two-particle Green functions. 
We found [18| that anticorrelated movements of particles 
govern the dynamics of particles inside Mott lobes, as the 
particles of different types exchange their positions. The 
two-particles anti-correlated Green function 

Ga(R) = ^ J2 ( a L a 0r+R a 0raAr+R + H ' C -) . ( 6 ) 
r 

measures the phase coherence of such exchange move- 
ments as a function of distance. Due to its defini- 
tion, G a cannot be larger that pop\ and is equal to 
G a (R) = Go(R)Ga(R) if there is no correlation between 
the movements of and A particles. 

In two dimensions, at low but finite temperatures, we 
expect, in some cases, to observe a Berezinsky-Kostcrlitz- 
Thoulcss (BKT) phase transition and the different Green 
functions to adopt a power law behavior at large distance 
R, G(R) oc R~ n , characteristic of the appearance of a 
quasi long range order (QLRO) in the phase, long range 
order (LRO) being achieved only at T = (rj varying be- 
tween 1/4 and as T is lowered) (26| . In other words, at 
finite T we can expect a superfluid where the phase is stiff 
but not ordered and not a BEC whith an ordered phase. 
At high temperature, the Green functions are of course 
expected to decay exponentially. On the finite size sys- 
tems which we study (L < 14), it is difficult to distinguish 
the QLRO from a true LRO, whereas one can easily dis- 
tinguish between the QLRO and exponentially decreas- 
ing regimes. This difficulty of distinguishing QLRO from 
LRO is also encountered in experiments where the sizes 
of systems that can be studied are typically of the same 
order as in our QMC simulations (hundreds of particles). 
The finite size results are, therefore, directly relevant to 
experiments. 

To elucidate the properties the model, we formulate 
it in terms of spins using a Schwinger bosons approach 
[27| . Defining the spin operators S* = (h 0r — fiAr) /2, 



= a 0r a\ r ^ and S r = a A r a 0r the Hamiltonian takes 
the form 

H = /U^E^r-l) + ^E«r (7) 
r r 

+ C/ 2 ^((5?) 2 -[(^) 2 + (^) 2 ]) (8) 

r 

where K is the hopping term Eq. ([T]). The terms in ([7]) 
are invariant under spin rotations. On the other hand, 
term Q favors pseudo-spin correlations to develop along 
the x axis if C/2 < or in the yz plane for U2 > 0. We 
note that the total spin S 2 = s r (s r + 1) on a given site 
is not fixed but depends on the total number of particles 
on the site (s r — n r /2) and will then fluctuate with this 
number. 

An order along z is measured through the densities 
or through density-density correlations of or A parti- 
cles, i.e. it corresponds to the polarization of the sys- 
tem. An order along the x or y axes is exposed through 
the behavior of G a which, in terms of spins, is equal to 
G a (R) = E r (^+R + Sy r S» +B )/L 2 . Our QMC algo- 
rithm allows the calculation of G a (R) but does not give 
access to correlations along the x and y axes indepen- 
dently 

We remark that the y and z axis have the same behav- 
ior, which means that we expect a spin QLRO in the yz 
plane to appear at low enough temperature for t/jj > 
in addition to the expected QLRO of the global phase 
of the particles discussed earlier. On our finite size sys- 
tems, this means that we should simultaneously observe 
a polarization and a QLRO for G a . We will call such 
a phenomenon a "quasi-polarization" (QP). For U<x > 
and finite T, we expect a QLRO to develop along the x 
axis and, consequently, no polarization but still a plateau 
at long distances in the function G a . 

III. U 2 > CASE 

At zero temperature, the results obtained with QMC 
and MFT were in good qualitative agreement [Tsj . The 
phase diagram, studied for densities up to two, exhibits 
three phases. The first two phases are incompressible 
Mott phases obtained for integer densities p = 1 and 
p = 2 for large enough interactions Uq. At zero tempera- 
ture, the entire p = 1 Mott phase is polarized, that is the 
system sustains a spontaneous symmetry breaking and 
the density of one type of particles becomes dominant. 
This polarization can be understood in the framework of 
an effective spin-1/2 model [lj| as the coupling in the 
yz plane is stronger than along the x axis. As expected, 
exchange movements of particles coexist with polariza- 
tion in this phase and G a develop a long range phase 
coherence. The p = 2 Mott phase is unpolarized and 
show no sign of exchange moves. In terms of spin, this 
corresponds to the fact that, neglecting the kinetic term, 
the ground state for a given site is uniquely determined 
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FIG. 1: (Color online) The limits of the p = 1 and p = 2 Mott 
regions (MI) for different values of the temperature kT and for 
different sizes L. As the temperature increases, the Mott lobes 
shrink and totally disappear for kT > 0.04 Uo- Outside the 
Mott, the system crosses over to a normal Bose liquid (NBL) 
and, eventually, transitions to a quasi-polarized superfluid 
(QP-SF) phase (see below). In the limit of zero temperature 
(here represented by the low temperature kT/Uo = 0.006), 
there is a direct phase transition between a Mott phase and 
a superfluid. 



as the state with = 0. The third possible phase is a 
polarized superfluid (SF) and occurs at any Uq when the 
density is incommensurate and also at small Uq for com- 
mensurate values. It is not possible to discuss this phase 
in terms of a simple effective spin degree of freedom as 
the number of particles on a site is not fixed. 

At zero temperature, the transition from the p = 1 
Mott phase to the SF is continuous, whereas at the tip 
of the p = 2 Mott lobe, the transition to the SF is first 
order for small values of U2/U0 becoming second order 
for larger values. This was predicted by the MFT (Tl| 
and confirmed by QMC simulations fill ]. 

A. Phase diagram at T / 

To map the phase diagram at finite temperature with 
QMC, we determine the limit of the Mott Insulator re- 
gions by measuring the density as a function of p and 
determining the boundaries of the plateaux indicating 
the almost incompressible regions. As explained in Sec. 
IIA, although the compressibility of these regions is very 
small, they are not strictly incompressible due to ther- 
mal fluctuations. The evolution from MI to NBL does 
not show any singularity and is then simply a crossover 
between two different limiting behaviors, incompressible 
and compressible, of the same unordered phase. Strictly 
speaking, a truly incompressible Mott phase exists only 
at T = but, following convention, we will continue to re- 
fer to this finite T region as a MI. To define the crossover 



boundary between the MI and the NBL shown in phase 
diagrams, we use the following criterion: when the den- 
sity deviates by 1% from the total integer density we 
consider that the system is no longer in the MI region 
but in the NBL region. 

Figure [1] shows the limits of the MI regions for differ- 
ent temperatures. As expected, the MI regions progres- 
sively disappear as the temperature is increased. The 
limits of the superfluid phase are located by direct mea- 
surement of p s . When the system is compressible and 
has zero superfluid density, there is a normal Bose liq- 
uid. Fig. [2] (top) shows the QMC phase diagram of the 
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FIG. 2: (Color online) Phase diagram of the system at 
kT = 0.03 Uo for U2/U0 = 0.1. The QMC simulations (top) 
show that a liquid (NBL) region appears between the quasi 
polarized superfluid (QP-SF) phase and the Mott region. The 
MFT (bottom) is not able to reproduce correctly this result 
as it does not take into account the kinetic term in the un- 
ordered (Mott or NBL) phase. It predicts first order phase 
transition between the SF and the other phase at the tip of 
the p = 1 and p = 2 lobes: the region of coexistence of the 
two phases is limited by the dashed line and the redline with 
squares. 
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FIG. 3: (Color online) The QMC phase diagram for p — 1 and 
for different sizes showing the quantum phase transition point 
between the Mott insulating phase (MI) and the superfluid at 
kT = 0. The finite T superfluid phase does not have BEC in 
this two-dimensional system. The system is quasi-polarized 
throughout the superfluid phase, but it is not in the NBL 
and MI. The Mott phase is polarized at T = and should be 
quasi-polarized at extremely small T. The transition between 
the SF and the NBL is continuous and there is only a crossover 
separating NBL from MI. 

system at a constant finite temperature. As the SF and 
MI regions are destroyed due to thermal fluctuations, an 
intermediate NBL region appears. The MFT used with 
success at T = 0, where it reproduces qualitatively the 
phase diagram, is unable to do so at finite temperature 
as explained in Sec. Ill Al Figure [2] (bottom) shows the 
MFT phase diagram where the afore mentioned problem 
clearly appears: the boundaries between the MI and NBL 
do not depend on the value of t/Uo and the MFT is un- 
able to give correct predictions regarding this crossover. 
However, the boundaries of the SF are reasonably well 
reproduced. A surprising result is that the transition be- 
tween the unordered phase and the SF appears to be first 
order at the tip of the p = 1 and p = 2 lobes in this MFT 
approach. At zero temperature, only p — 2 showed a 
first order transition. As will be shown below, for finite 
temperatures, this is in total contradiction with the re- 
sults obtained by QMC simulations that show continuous 
phase transitions for p = 1 and p = 2. So this MFT pro- 
vides incorrect description of the phase transition and of 
the position of the different regions at finite temperature. 

For fixed integer density at zero temperature, we have 
a quantum phase transiti on ( QPT) between the MI and 
SF phases. As expected [28j . when the temperature is 
increased from zero, an intermediate compressible un- 
ordered region appears between the superfluid phase and 
the Mott region, namely the normal Bose liquid region 
(NBL). This is observed for p = 1 (Fig. |3J) and p = 2 
(Fig.©. 

As for the possible polarization of the different phases, 
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FIG. 4: (Color online) The QMC phase diagram for p — 2 
and for different sizes showing the quantum phase transition 
point between the MI and the SF occurring only at kT — 0. 
The superfluid phase is always quasi-polarized whereas the 
NBL and MI are not. The SF-NBL transition is continuous, 
whereas at T = the MI-SF quantum phase transition was 
found to be discontinuous for U2/U0 small enough The 
dash line corresponds to the case studied in Fig. El 

we observe that the SF phase appears to be polarized at 
finite T as it is at T = 0: the histogram of the density 
of one of the species shows two peaks at low and high 
densities (see Fig. [5]). However, due to the continuous 
symmetry in the yz plane of the pseudo-spin part of the 
Hamiltonian Eq. (|5J|, no LRO exists at finite T. There- 
fore, this apparent polarization is due to the finite size 
of the system, and is, in fact, a quasi-polarization. An 
intuitive way to understand this quasi-polarization is as 
follows: In the superfluid phase, the pseudo-spins in the 
yz plane are stiff and appear to be mostly aligned in the 
same direction on a finite lattice, such as the case here. 
But since the symmetry is not broken, this "magnetiza- 
tion" direction in the yz plane will drift and point in all 
directions. As the direction changes, so does the projec- 
tion of this pseudo-spin on the z-axis. In other words, the 
polarization drifts too, and changes with time, giving the 
double peak structure to the polarization histogram, Fig. 
[5] We note that experimental systems are typically of the 
same sizes as the ones we study here and, therefore, this 
polarization drift will be present in these experiments 
too: The particle content of the system will change as a 
function of time. This quasi-polarization disappears as 
T increases when the system undergoes a thermal BKT 
phase transition into the NBL. In other words, the entire 
SF phase is quasi-polarized but the NBL is not (see Fig. 
ED 

A histogram of the polarization in the p = 1 MI region 
shows that as soon as the temperature is increased from 
zero, the polarization disappears and the populations be- 
come balanced. According to the effective pseudo-spin 
model [H, , it is possible to observe quasi-polarization 
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FIG. 5: (Color online) Probability P(pa) as a function of 
temperature for L = 8, U 2 /U = 0.1, p = 2 and t/U = 0.08. 
When kT/Uo < 0.175, the system is in the quasi-polarized SF 
phase. For kT/Uo > 0.175, the system is in the unordered 
unpolarized phase. The temperature at which the quasi- 
polarization disappears is approximately the same as the tem- 
perature where p s becomes zero (see the vertical dashed line 
in Fig. UJ. Only the SF phase is quasi-polarized. 



as in the SF phase. In this case, the coupling generating 
these spin correlations and the polarization of the MI is 
of order t 2 /Uq. However, even for temperatures as low as 
kT/Uo — 0.01, the system is already in the regime where 
correlations decay exponentially and we do not observe 
any QLRO at finite T in the Mott region. This is con- 
firmed by the behavior of the Green functions that is 
shown in Fig. [S] In Fig. GHa) we show the anticorrelated 
Green function G a , Eq. , in the p = 1 MI region as T 
is increased. We see that as soon as T becomes finite, G a 
decays exponentially with distance, contrary to its con- 
stant value at long distance observed at T = . This 
exponential decay of course persists in the NBL region. 

In Fig. [Bfb) we show G A (R), Eq. ©, for various T 
values at p = 2 and t/Uo = 0.08 which, at T = 0, puts the 
system in the SF phase. We see that for low T, G\ decays 
slowly. This decay is expected to be a power law but the 
system size is too small to show that unambiguously. As 
T is increased, the system transitions into the NBL where 
Ga exhibits exponential decay clearly distinguishing the 
NBL and SF phases. A similar behavior is found for Go 
and G a , the latter being expected since it accompanies 
the presence of quasi-polarization. 



B. Nature of the transitions 

At finite temperature, kT > 0.005 Uq, the transitions 
between the SF and the NBL are continuous for p = 1 as 
well as for p = 2 (and of course for all other densities). 
Since this transition takes place at constant density and 
is therefore a phase-only transition, it is expected to be 
in the Berczinsky-Kostcrlitz-Thouless universality class 
(BKT) for our two-dimensional system. Actually, below 
the critical temperature, we have two different quasi long 



10" 



10 



10' 



10- 



10 



10 



(a) p=1,t/U.=0.02 



g„(R) 






kTAJ 


-0.0004 


e-o 


kT/U 


-0.0100 


O-D 


kT/U 


-0.0200 


0-0 


kT/U 


-0.0300 


G-O 


kT/U 


-0.0400 


e-H 


kT/U 


-0.0500 




kT/U 


-0.0700 




kT/U 


-0.1000 


B-H 


kT/U 


-0.1200 


<^kT/U 


-0.1400 


G-O 


kT/U 


-0.1600 


O-H 


kT/U 


-0.1800 


G-e 


kT/U 


-0.2000 




FIG. 6: (Color online) Evolution of the Green functions G a , 
Eq. ©, in the p = 1 MI and Ga, Eq. ©, in the SF as func- 
tions of distance R for different temperatures, with L = 12 
and U2/U0 = 0.1. (a) In the p = 1 MI, G a decays exponen- 
tially at finite T, unlike its T = behavior where it reaches a 
constant value [l8|]. A regime of QLRO is not observed and 
would occur only at very small T. The exponential decay 
persists in the NBL at higher T. (b) Ga at p = 2 for T values 
taking the system from the SF to the NBL. Same case as in 
Fig. O In the SF, Ga decays slowly as a power law, but de- 
cays exponentially in the NBL. Similar behavior is found for 
Go and G a . 



range orders that occur: the global U(l) phase QLRO 
associated with the superfluid behavior and the pseudo- 
spin one, associated with the quasi-ordering of the spins 
in the yz plane, i.e. the so-called quasi polarization. As 
explained previously, these two QLRO appear simultane- 
ously. Since the transition is BKT, we first determined 
the transition temperature using the universal jump of 
the superfluid density (29| . where, at the transition tem- 
perature, T c , we have p s {T c ) = kT c /irt. To observe this, 
we calculated p s as a function of temperature and deter- 
mined T c (p s ) graphically as the intersection of p s (T) with 
kT/irt (see Fig. [7]). We also calculated the specific heat 
G and determined the transition temperatures T c (C max ) 
as the temperature where C reaches its maximum (Fig. 
[7|. Our QMC simulations were done for L = 8 because 
G is extremely difficult to calculate at low temperature 
for U2 > for larger sizes. 

Tabic |T] compares the values of T c obtained from the 
universal jump and from the maximum of G. The values 
are in agreement confirming the universal jump hypothe- 
sis and the BKT nature of the transition and determining 
the transition temperature for the studied size. We did 
simulations for sizes up to L = 14 to examine the effect 
of finite size for this transition (see Fig. As expected, 
the transition gets sharper with increasing size. However, 
it is too difficult to obtain results with small enough error 
bars for G on these large sizes. 

As mentioned earlier, the evolution from the MI region 
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Filled symbols: p s 
Open symbols: C/C 



A^t/U =0.06 

■ — ■ t/U =0.07 




Grand canonical ensemble 




FIG. 7: (Color online) Superfluid density p s and specific heat 
C as functions of temperature for different values of t/Uo. 
The universal jump condition states that p a (T c ) = kT c /irt. 
The transition temperatures are calculated independently by 
determining the maximum of C(T). 



t/Up T c (£s) T c (Cmax) 

0.06 1.49 ±0.10 1.67 ±0.15 

0.07 1.82 ±0.05 1.85 ±0.15 

0.08 1.99 ±0.05 2.00 ±0.15 

0.09 2.10 ±0.05 2.10 ±0.20 

0.10 2.19 ±0.05 2.20 ±0.20 

TABLE I: Values of the transition temperatures for different 
values of t/Uo for an 8 x 8 system. T c (p 3 ) is obtained from the 
universal jump while T c (C ma x) from the maximum of C(T). 
The values are in agreement and show that the universal pre- 
diction is verified. The data are taken from Fig. [7] 



Transition SF-NBL 

AF, U 2 =0.1U , p=2 



An*L=8, t/U =0.06 
«_OL=10, t/U =0.06 
OOL=12, t/U =0.06 
rL=14, t/U =0.06 
L=8, t/U =0.10 
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O-oL=12,t/U =0.10 
v^L=14, t/U =0.10 




FIG. 8: (Color online) Superfluid density p s as a function of 
T for t/Uo = 0.06 and 0.10 and for different sizes L. The 
transition becomes sharper as L is varied from L = 8 to L — 
14. 



FIG. 9: (Color online) p, p B and k as functions of p. As one 
goes from the MI to the NBL, the density varies continuously 
indicating a crossover between the two regions. The two re- 
gions are distinguished by the integer/non-integer value of the 
density and by the larger value of the local compressibility k 
in the liquid. The superfluid phase has nonzero p s and a much 
larger local compressibility than both the MI and NBL. k has 
been multiplied by an arbitrary factor A for better visibility. 



to the NBL is a continuous crossover. A plot of the den- 
sity as a function of p shows no sign of a first order phase 
transition in the form of a jump in the density as one ap- 
proaches the Mott plateaux (see Fig. [9J. There is no 
phase transition between MI and NBL since, in addition 
to the absence of a first order transition, no symmetries 
are broken. Then, at finite temperature, there is only a 
crossover between the MI and the NBL and not the phase 
transition predicted by MFT. 

At zero temperature, we have shown [l8| that the 
Mott-SF transition is always second order for p = 1 but 
is first order near the tip of the p — 2 Mott lobe when 
U2/U0 is small enough, for example U2/U0 = 0.1. Hence, 
while it is easy to imagine that the continuous NBL-SF 
transitions observed at moderate temperatures persists 
at low temperature for p = 1, the p = 2 case where the 
behavior is different at zero and finite temperatures re- 
quires a separate study for the low temperature regime. 
Probing temperatures as low as kT = 0.005[/n (see Fig. 
1101) . the transition still appears continuous, which is to be 
compared with the temperature at which the MI region 
is destroyed kT « 0.05t/ (see Figs. |3]and|4]). Although 
it cannot be excluded that the transition is discontinuous 
in a small range of temperatures, this range would be ex- 
tremely narrow. In order to observe the first order QPT 
in our previous work [Lsl ] we used temperatures that were 
of order kT/Uo ~ 10~ 3 . The transition then appears dis- 
continuous only at extremely low temperatures. 
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FIG. 10: (Color online) Transition from the normal Bose liq- 
uid to the superfluid phase at low temperatures. The tran- 
sition always appears continuous for low temperatures. The 
discontinuity is only found in the zero temperature regime 
(kT < W- 3 U ). 



IV. U 2 <0 CASE 

It was shown for U2 < [TH that all the phases are 
unpolarized and that transitions between the different 
phases are all continuous in the T = limit. The coher- 
ent anticorrelatcd movements are present in the p = 1 
and in the p = 2 Mott phases. 

Proceeding as in the U2 > case, we determine the 
boundaries of the Mott regions at finite T, Fig. [TTJ by 
calculating the density and the boundary of the super- 
fluid region by measuring p s . In this case, we chose not 
to present results from MFT since it shows the same lim- 
itations as in the U2 > case. We obtain the phase 
diagram for p = 1 (shown in Fig. I12[) and a similar one 
for p = 2 (not shown here) . Similarly, we used histograms 
of the density similar to Fig. [5] to confirm that all these 
phases remain unpolarized at finite temperature as they 
are at zero temperature. The absence of polarization at 
T can be understood qualitatively from Eq. (|5J). For the 
present case, U2 < 0, the last term in Eq. (J5J favors 
the alignment of the spins along the x-axis in the low 
T phase. Consequently, the polarization, S z . is always 
zero. 

As in the U2 > case we observe a slow decay of the 
Green functions G a in the superfluid phase and the tran- 
sition is shown to be of the BKT type using the universal 
jump argument (see Table HI| . 

We also examined the anticorrelated Green function 
G a in both p = 1 and p = 2 MI regions at various temper- 
atures. In the zero temperature limit, both these phases 
exhibit nonzero values of G a at long distances. This is ex- 
pected in both cases, considering the pseudo-spin Hamil- 
tonian. Neglecting the kinetic energy, for U2 < 0, the en- 
ergy is minimized by maximizing (S^) 2 on each site. For 



t/Uo T c (p s ) T c (C m ax) 

0.06 1.85 ±0.05 1.75 ±0.10 

0.07 2.00 ±0.05 2.00 ±0.10 

0.08 2.10 ±0.05 2.05 ±0.05 

0.09 2.14 ±0.05 2.16 ±0.05 

0.10 2.18 ±0.05 2.15 ±0.05 

TABLE II: Values of the SF to NBL transition temperatures 
for different values of t/Uo and for a L — 12 size for U2/U0 = 
-0.1 and (0 = 2. Similar to Table [fl 
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FIG. 11: (Color online) The limit of the p = 1 and p = 2 
Mott regions (MI) for different values of the temperature kT 
and for different sizes L in the U2 < case. As in the positive 
U2 case, the Mott phases disappear for kT > 0.04J7o- In the 
T — limit, there is a direct transition between the MI and 
SF. 



0.2 



0.15 - 



h- 



0.1 - 



0.05 



A L=8 
□ L=10 
00 |_=12 



U 2 /U =-0.1,p=1 



NBL 



Ml 




Superfluid 
unpolar 



0.05 



t/lL 



0.1 



0.15 



FIG. 12: (Color online) The phase diagram for p = 1 in the 
U2 < case. All phases are unpolarized and the SF-NBL 
transition at finite temperature as well as the MI-SF tran- 
sition at zero temperature are continuous. A similar phase 
diagram is observed for p = 2. 
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FIG. 13: (Color online) The anticorrelated Green functions 
in the p = 1 (a) and p — 2 (b) Mott regions for negative 
Ui. In both cases, the long range order of the anticorrelated 
movements that was present at T = is lost even for very 
small temperatures. 



p = 1 and p = 2, there are two degenerates states that 
achieve this: the S* = ±1/2 for p = 1 and the S* = ±1 
for p = 2. These degenerate ground states are coupled by 
second order hopping processes which lift this degener- 
acy through coherent anticorrelated movements [18j. In 
Fig. [T31 we show G a in the p = 1 and p — 2 MI regions 
for different temperatures. As expected, the phase co- 
herence once again completely disappears rapidly as the 
temperature is increased from zero and we do not observe 
any sign of quasi long range phase coherence in the MI 
region, as well as in the NBL phase. 



V. CONCLUSION 

Studying a bosonic spin-1/2 Hubbard model at finite 
temperature and comparing to the T = case, we find 
that the effect of temperature is dramatically different de- 
pending on the phase we consider. The superfluid phases 
are essentially unchanged by raising the temperature: 
The long range order present at T = is transformed 
into QLRO. On the other hand, the MI regions are dras- 
tically modified as the polarization that occurs in certain 
cases is almost immediately wiped out by thermal fluc- 
tuations: low temperature quasi polarized states are not 
found in the regime of temperatures we studied. 

This is due to the fact that the energy scales associated 
with the polarization of the system take different values 
in these different phases. In the MI regions, the polariza- 
tion is due to the coupling between different low energy 
degenerate Mott states, as emphasized by the pseudo- 
spin theory [lj|. These couplings arc of order t 2 /Uo and 
the pseudo-spin quasi-ordering vanishes very rapidly with 



temperature. On the other hand, in the superfluid, the 
energy scale associated with the coupling of pseudo-spins 
is obviously much larger. While it is not possible to spec- 
ify this scale as precisely as in the Mott phases, due to the 
itinerant nature of the particles in the superfluid regime, 
a simple argument shows that it is of order L^: when- 
ever the particles enter the superfluid phase and adopt 
delocalized states, they overlap; there is an interaction 
cost which is then of order Uo for identical particles and 
Uq + U2 for different ones. This favors having the par- 
ticles mostly of the same type for U2 > and leads to 
a quasi polarization. For C/ 2 < 0, the interaction favors 
having a mixture of particles and gives a system without 
any sign of polarization. In both cases, the energy scale 
is typically | C/2 1 

The finite temperature phase diagram presented here 
is important for the proper interpretation of experimen- 
tal realization of this and related systems using ultra- 
cold atoms loaded on optical lattices. Such experiments 
are, of course, always at finite temperature. Similarly 
to what happens for fermions, we observe that the spin 
correlations in the Mott phases will be very difficult to 
access experimentally as the associated energy scales are 
very small and as the correlations are almost immediately 
wiped out by thermal fluctuations. On the other hand, 
due to the relatively small sizes of experimental systems 
and to its larger associated energies, the quasi-ordering of 
spins in the superfluid phase should be immediately visi- 
ble experimentally and indistinguishable from true polar- 
ization of the system. In a finite size system, one would 
expect to observe a slow drift of the polarization as the 
yz symmetry of the system is restored over time. 

From a more technical point of view, we have eluci- 
dated the limitations of the MFT commonly used in the 
literature. We found that it is unable to distinguish cor- 
rectly the MI and NBL regions. The MFT does predict 
reasonably well the NBL-SF boundaries but not the na- 
ture of the transition which is sometimes predicted to be 
of first order whereas direct QMC simulations, and sym- 
metry considerations, show that it is in the BKT univer- 
sality class. Furthermore, we found that MFT predicts 
direct first and second order transitions at finite T be- 
tween the MI and SF phases which QMC shows do not 
exist, (see Figs. EJHH2])- Similar incorrect MFT behav- 
ior was found for the spin-1 model and the same caveats 
should be applied for example in Ref . [2fj| . 



Acknowledgments 

This work was supported by: the CNRS-UC Davis 
EPOCAL LIA joint research grant; by NSF grant OISE- 
0952300; an ARO Award W911NF0710576 with funds 
from the DARPA OLE Program. We would like to thank 
Michael Foss-Feig and Ana-Maria Rcy for useful input 
and discussion. 



10 



[1] R.P. Feynman, Int. J. Theor. Phys. 21, 467 (1982). [15 
[2] I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. [16 

80, 885 (2008). [17 
[3] M. Greiner, O. Mandel, T. Esslinger, T.W. Hansen, and 

I. Bloch, Nature 415, 39 (2002) 
[4] G.-B. Jo, Y.-R. Lee, J.-H. Choi, C.A. Christensen, T.H. [18 

Kim, J.H. Thywissen, D.E. Pritchard, and W. Ketterle, 

Science 325, 1521 (2009). 
[5] A. Griesmaier, J. Werner, S. Hensler, J. Stuhler, and T. [19 

Pfau, Phys. Rev. Lett. 94, 160401 (2005). 
[6] M.W. Zwierlein, A. Schirotzek, C.H. Schunck, and W, [20 

Ketterle, Science 311, 492 (2006); 
[7] Y. I. Shin, A. Schirotzek, C.H. Schunck, and W. Ketterle, [21 

Phys. Rev. Lett. 101, 070404 (2008). 
[8] M. Vengalattore, S.R. Leslie, J. Guzman, D.M. Stamper- [22 

Kurn, Phys. Rev. Lett. 100, 170403 (2008); 

M. Vengalattore, J. Guzman, S. R. Leslie, F. Ser- [23 

wane, and D.M. Stamper-Kurn, Phys. Rev. A 81, 053612 

(2010). [24 
[9] D.M. Stamper-Kurn and W. Ketterle, in Coherent 

Atomic Matter Waves, edited by R. Kaiser, C. West- [25 

brook, and F. David, Springer, p. 137 (2001). 
[10] G.G. Batrouni, V.G. Rousseau, and R.T. Scalettar, [26 

Phys. Rev. Lett. 102, 140402 (2009). 
[11] K.V. Krutitsky and R. Graham, Phys. Rev. A 70, 063610 [27 

(2004); K.V. Krutitsky, M. Timmer and R. Graham, 

Phys. Rev. A71, 033623 (2005). 
[12] J. Larson and J.-P. Martikainen, Phys. Rev. A80, 033605 [28 

(2009). 

[13] W. V. Liu and C. Wu, Phys. Rev. A74,13607 (2006). [29 
[14] C. Wu, W. V. Liu, J. E. Moore and S. Das Sarma, Phys. 
Rev. Lett. 97 190406 (2006). 



C. Wu, Modern Physics Lett ers B23 1 (2009) . 
M. Foss-Feig and A.-M. Rev. larXiv:1103.0245k 2. 

L. de Forges de Parny, M. Traynard, F. Hebert, V.G. 
Rousseau, R.T. Scalettar, and G.G. Batrouni, Phys. Rev. 
A 82, 063602 (2010). 

L. de Forges de Parny, F. Hebert, V.G. Rousseau, R.T. 
Scalettar, and G.G. Batrouni, Phys. Rev. B 84, 064529 
(2011). 

A. B. Kuklov and B. V. Svistunov, Phys. Rev. Lett. 90, 
100401 (2003). 

R. V. Pai, K. Sheshadri, and R. Pandit, Phys. Rev. B 
77, 014503 (2008). 

T. Kimura, S. Tsuchiya, and S. Kurihara, 
Phys. Rev. Lett. 94, 110403 (2005). 

S. Takayoshi, M. Sato, and S. Furukawa, Phys. Rev. A 
81, 053606 (2010). 

V.G. Rousseau, Phys. Rev. E 77, 056705 (2008) ; 
V.G. Rousseau, Phys. Rev. E 78, 056707 (2008). 

D. M. Ceperley and E.L. Pollock, Phys. Rev. B39, 2084 
(1989). 

M. Eckholt and T. Roscilde, Phys. Rev. Lett. 105, 
199603 (2010). 

M. Le Bellac, "Quantum and Statistical Field Theory", 
Oxford University Press (1992). 

A. Auerbach and D.P. Arovas in "Introduction to Frus- 
trated Magnetism" , edited by C. Lacroix, P. Mendels, and 
F. Mila, Springer (2011). 

S. Sachdev, "Quantum Phase Transitions", Cambridge 
University Press (1999). 

D.R. Nelson and J.M. Kosterlitz, Phys. Rev. Lett. 39, 
1201 (1977). 



